COMMENT ⊗   VALID 00007 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	begin "AWFIL"
C00003 00003	! requires and other declarations
C00005 00004	! local procedures
C00008 00005	! routines to fetch input
C00010 00006	! main program here
C00011 00007	end "AWFIL"
C00012 ENDMK
C⊗;
begin "AWFIL"

require "HEADER.SAI[LIB,ROB]" source_file;

! lotsa stuff drawn from
	JAMLIB.JAM[UP,DOC]
	JAMLIB.SAI[SUB,SYS]
	JAMLIB.HDR[SUB,SYS];

! This program takes the frequency response of a system and calculates the
  coefficients for a set of second order filters which would flatten the
  response of the system.

  The input parameters are:
	something that describes the magnitude versus frequency (T.B.A.)
	the number of second order sections available

  The outputs are:
	various plots of the system response
	the coefficients of the second order sections

  ;
! requires and other declarations;

require "JAMLIB.REL[SUB,SYS]" library;
    external integer procedure POLY(
	real array A; reference real array ZR,ZI; integer N);
    external procedure RTRAN(
	real array X,A,B; integer log2N; boolean inverse);
    external real procedure AUTSOL(
	real array R,A,K; integer M);
    external string procedure DEFAULT(string Inp; reference string Dev;
	string Devext,Devdev);
    external boolean procedure SEGSYN(string dev,file;
	reference record_pointer(any_class) L; reference integer nfns);

    record_class seg(integer type; record_pointer(seg) next,last;
	string name;
	integer nsegs;
	real mintime,maxtime,minval,maxval;
	real array times,values);


require "MAGFRG.TXT" SOURCE_FILE;	! To preload MagFrq array;
real array MagFrq[0:NPoints];

! local procedures;

procedure COMBINE(
    real array Rls,Ims,Ai,Bi; integer M);
    begin "COMBINE"
	boolean array taken[1:M];
	integer ind,upto,jj;

	arrclr(taken);
	ind←1;
	for upto←1 step 1 until M do
	begin "FFCR"
	  if taken[upto] then continue "FFCR";
	  if abs(Ims[upto])<.0000001 then
	  begin "ISR"
	    for jj←upto+1 step 1 until M do
	      if abs(Ims[jj])<.0000001 then done;
	    if jj≤M then
	    begin "ISANO"
		Ai[ind]←-Rls[upto]-Rls[jj];
		Bi[ind]←Rls[upto]*Rls[jj];
		taken[jj]←true;
	    end "ISANO"
	    else begin "NOTAN"
		Ai[ind]←-Rls[upto];
		Bi[ind]←0;
	    end "NOTAN";
	  end "ISR"
	  else begin "ISI"
	    for jj←upto+1 step 1 until M do
	      if abs(Ims[upto]+Ims[jj])<.0000001
		∧ abs(Rls[upto]-Rls[jj])<.0000001 then done;
	    if jj>M then usererr(0,0,"Complex filter?????");
	    Ai[ind]←-2*Rls[upto];
	    Bi[ind]←Rls[upto]↑2+Ims[upto]↑2;
	    taken[jj]←true;
	  end "ISI";
	  ind←ind+1;
	end "FFCR";
    end "COMBINE";

procedure ArrSqr(real array Arr);
  ⊂ "ArrSqr"
  define ndims = {-1}, nwrds = {0}, arlo = {1}, arhi = {2};
  integer LoBound, HiBound, I;
  if arrinfo(Arr,ndims) = 1
    then
      ⊂
      LoBound ← arrinfo(Arr,arlo);
      HiBound ← arrinfo(Arr,arhi);
      for I ← LoBound step 1 until HiBound do Arr[I] ← Arr[I]↑2;
      ⊃
    else
      outstr("ArrSqr: multiply dimensioned array.  No good.  No change."&↓);
  ⊃ "ArrSqr";

! routines to fetch input;

record_pointer(seg) procedure GetFns;
  ⊂ "GetFns"
  while true do
    ⊂ "floop"
    define DefDev = {"DSK"}, DefExt = {"FUN"};
    string FilNam,DevNam;
    record_pointer(seg) FunLst, lp;
    integer NFuncs;

    outstr("."&DefExt&" file for frequency response: ");
    FilNam ← DEFAULT(inchwl,DevNam,DefExt,DefDev);
    if SEGSYN(DevNam,FilNam,FunLst,NFuncs) then done;
    outstr("Non-existant or screwy function file (?), try again..."&↓);
    ⊃ "floop"
  outstr(CVS(NFuncs)&" functions found:"&↓);
  lp←Funlst;
  while lp do
    ⊂
    outstr(seg:name[lp]&TAB);
    lp←seg:next[lp];
    ⊃;
  outstr(↓);
  return(FunLst);
  ⊃ "GetFns";

record_pointer(seg) procedure GetSeg(record_pointer(seg) LP; string NAME);
  ⊂ "GetSeg"
  ! searches down the list until it finds a SEG function named NAME, and
    returns a pointer to it.  Returns null_record if not found;
  while LP do
    ⊂
    if equ(NAME,seg:name[LP]) then done;
    LP ← seg:next[LP];
    ⊃;
  return(LP);

! main program here;

⊂ "MAIN"
! square the magnitude response;
ArrSqr(MagFrq);

! take the inverse FFT of the squared magnitude response (yeilds autocorrelation);

! run a Linear Predictor over it to get the filter coefficients;

! factor the filter coefficients into their component roots;

! combine the roots into conjugate pairs;

! output the roots, now suitable for a cascaded second order system;

end "AWFIL";